#include <math.h>
#include <vector>
#include <stdlib.h>
#include <iostream>
#include <limits>
#include <cstdio>
#include <cstdlib>
#include <string.h>
#include <float.h>
#include <fstream>
#include "geometry.h"
#include <boost/numeric/interval.hpp>


double EPSILON = std::numeric_limits<double>::epsilon();

////////// Simple //////////
int turn(point a, point b, point c)
{
   double res = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
   double modres = abs((b.x - a.x) * (c.y - a.y)) + abs((b.y - a.y) * (c.x - a.x));
   double check = 4 * EPSILON * modres;
   if(res > check)
      return 1;
   if(res < -check)
      return -1;
   return 0;
}
////////// Simple //////////


///////// Interval //////////
int interval_turn(point a, point b, point c)
{
   using namespace boost::numeric;
   using namespace interval_lib;

   interval<double> iax(a.x), iay(a.y),
                    ibx(b.x), iby(b.y),
                    icx(c.x), icy(c.y);

   interval<double> res = (ibx - iax) * (icy - iay) - (iby - iay) * (icx - iax);
   if(res.lower() > 0)
      return 1;
   if(res.upper() < 0)
      return -1;
   return 0;
}
///////// Interval //////////


///////// Adaptive precision floating point arithmetic //////////
double sum(double a, double b, double & roundoff)
{
   double res = a + b;
   double bv = res - a;
   double av = res - bv;
   double br = b - bv;
   double ar = a - av;
   roundoff = ar + br;
   return res;
}

template <size_t s>
void split(double a, double & hi, double & lo)
{
   double c = ((1LL << s) + 1LL) * a;
   double ab = c - a;
   hi = c - ab;
   lo = a - hi;
}

double mul(double a, double b, double & roundoff)
{
   double res = a * b;
   double a_hi, a_lo, b_hi, b_lo;
   split<std::numeric_limits<double>::digits / 2 + std::numeric_limits<double>::digits % 2>(a, a_hi, a_lo);
   split<std::numeric_limits<double>::digits / 2 + std::numeric_limits<double>::digits % 2>(b, b_hi, b_lo);
   double e1 = res - (a_hi * b_hi);
   double e2 = e1 - (a_lo * b_hi);
   double e3 = e2 - (b_lo * a_hi);
   roundoff = (a_lo * b_lo) - e3;
   return res;
}

template <size_t m>
void grow_expansion(double * e, double b, double * h)
{
   double q = b;
   for(size_t i = 0; i < m; i++)
      q = sum(e[i], q, h[i]);
   h[m] = q;
}

template <size_t m, size_t n>
void expansion_sum(double * e, double * f)
{
   for(size_t i = 0; i < n; i++)
      grow_expansion<m>(e + i, f[i], e + i);
}

template <size_t n>
int ap_sign(double * e)
{
   for(int i = n - 1; i >= 0; i--)
   {
      if(e[i] > 0)
         return 1;
      if(e[i] < 0)
         return -1;
   }
   return 0;
}

int ap_turn(point a, point b, point c)
{
   double mult[12];

   mult[0] = mul(b.x, c.y, mult[1]);
   mult[2] = mul(-b.x, a.y, mult[3]);
   mult[4] = mul(-a.x, c.y, mult[5]);
   mult[6] = mul(-b.y, c.x, mult[7]);
   mult[8] = mul(b.y, a.x, mult[9]);
   mult[10] = mul(a.y, c.x, mult[11]);

   expansion_sum<2, 2>(mult, mult + 2);
   expansion_sum<2, 2>(mult + 4, mult + 6);
   expansion_sum<2, 2>(mult + 8, mult + 10);
   expansion_sum<4, 4>(mult, mult + 4);
   expansion_sum<8, 4>(mult, mult + 8);

   int res = ap_sign<12>(mult);

   return res;
}
///////// Adaptive precision floating point arithmetic //////////


//////// left turn predicate //////////
int left_turn(point a, point b, point c)
{
        int res = turn(a, b, c);
        if(res != 0)
                return res;

        res = interval_turn(a, b, c);
        if(res != 0)
                return res;

        res = ap_turn(a, b, c);
        return res;
}
//////// left turn predicate //////////
